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Abstract. We present a description of our numerical code BLAZAR. This code calculates spectra and light 
curves of blazars during outbursts. The code is based on a model in which the non-thermal flares in blazars 
are produced in thin shells propagating down a conical jet with relativistic velocities. Such shells may represent 
layers of a shocked plasma, enclosed between the forward and reverse fronts of an internal shock. In the model 
adopted by us, the production of non-thermal radiation is assumed to be dominated by electrons and positrons 
which are accelerated directly, rather then injected by pair cascades. The code includes synchrotron emission 
and inverse-Compton process as the radiation mechanisms. Both synchrotron photons and external photons are 
included as the seed photons for Comptonization. At the present stage, the code is limited to treat the inverse 
Compton process only within the Thomson limit and is specialized to model radiation production in the flat 
spectrum radio quasars. As an example, we present the results of modeling an outburst in 3C 279 - the most 
extensively monitored 7-ray - bright quasar. 
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1. Introduction 

One of the greatest achievements of the recently retired 
Compton Gamma-Ray Observatory (CGRO) is the dis- 
covery of gamma ray emission from a subclass of Active 
Galactic Nuclei (AGNs) known as blazars. More than 50 
such sources were detected at energies above lOOMeV 
(Mukherjee et al. I1997[l . Energy constraints and 7-ray 
absorption arguments require that radiation emission in 
blazars must be beamed (Mattox et al. Q.993). Indeed, 
it is now widely believed that the entire electromagnetic 
spectrum observed in these objects is dominated by non- 
thermal radiation produced in a jet pointing close to the 
line of sight (Dondi & Ghisellini ll995|l . Although roughly 
divided into two classes: flat spectrum radio quasars 
(FSRQ) and BL Lacs objects, all blazars share common 
characteristics: large-amplitude, rapid variability; smooth 
continuum emission in all observable bands; and high lin- 
ear polarization. 

Spectrum of a blazar consists of two broad components 
(von Montigny et al. 1995 ). The low energy component has 
a peak within IR-to- X-ray range and is usually attributed 
to Doppler-boosted synchrotron radiation. The high en- 
ergy component peaks in the MeV-TeV energy range and 
is very likely produced by inverse-Compton process. Both 
components are highly variable, with time scales rang- 
ing from years to a fraction of a day. Analysis of 7-ray 
light curves seems to suggest that variability patterns of 
blazars are a superposition of short term flares (Magdziarz 



et al. H997f) . The flares detected in different spectral bands 
appear to be correlated (Macomb et al. 119951 Wagner et 
al-HnSO Wehrle et al. UMfy . This behavior suggests co- 
spatial production of the high and low energy components 
and indicates that significant fraction of the jet energy is 
dissipated in the localized events at sub-parsec distances 
from the center. Such events are likely to arise in inter- 
nal shocks (Sikora et al. 1994; Spada et al. I2001[> or at 
the sites of the magnetic field reconnection (Romanova & 
Lovelace ITM21 Blackmail ESS}. 

Both the shocks and the reconnection sites provide fa- 
vorable conditions for efficient acceleration of particles. 
Unfortunately, the present theories of particle accelera- 
tion are still not sufficiently developed to provide quanti- 
tative predictions regarding such issues as what fraction 
of dissipated energy is used to accelerate particles, how 
much power is channeled into relativistic protons vs. rel- 
ativistic electrons/positrons, and what maximum energy 
can a particle gain. In particular, it is not possible to de- 
duce solely on theoretical grounds whether protons are 
accelerated to energies sufficient to support - via inelas- 
tic collisions with photons - the pair cascades, as is sug- 
gested by hadronic models (Mannhein & Biermann 1992), 
or if they are "radiatively" passive returning all energy 
gained during the dissipative events back to the flow via 
adiabatic expansion. In the latter case, represented by the 
leptonic models, radiation production in blazars is totally 
dominated by electrons and positrons accelerated directly. 
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There are several electron/positron acceleration mecha- 
nisms which are likely to operate in shocks (Levinson) 1996 
Hoshino et al. McClements et al. 113571 Shimada & 

Hoshino 2000 ; Hoshino & Shimada 2002 ) and reconnection 
sites (Larrabee et al. 120021 Lvutikov l2002fl . but such the- 
oretical considerations cannot as yet predict the spectra 
of accelerated particles and their minimum and maximum 
energies. 

Theories of radiation mechanisms are much better de- 
veloped. They are quantitative and for a given energy dis- 
tribution of the accelerated/injected particles, geometry 
and kinematics of the source, and external radiation held 
environment, one can make specific predictions regarding 
the radiation spectra. By confronting those theories with 
observations, it is possible to verify the model and deter- 
mine its parameters. Using this approach, one can already 
exclude some radiation scenarios. This concerns e.g. mod- 
els with pair cascades, which predict X-ray spectra that 
would be too soft as compared with observations of some 
FSRQ (Sikora & Madejski EMU)- Therefore, for these ob- 
jects, direct electron/positron acceleration mechanism is 
favored. 

In leptonic models, both the high-energy and the low- 
energy components of the spectrum are produced by the 
same population of relativistic electrons. At low ener- 
gies, the emission comes from the synchrotron process, 
while at higher energies, it is dominated by the inverse 
Compton scattering. There is still some debate about the 
source of seed photons for the inverse Compton process. 
The most obvious choice involves synchrotron photons 
from the low energy component. Models based on this 
assumption are called Self-Synchrotron-Compton (SSC) 
models. Originally proposed by Konigl ( 198TJ to explain 
production of X-rays and 7-rays (already detected by 
then from 3C 273 by COS-B [Swanenburg et a1. [T37B| ). 
currently the SSC models prove to be successful in ex- 
plaining general features of BL Lac objects (Ghisellini & 
Maraschi IT§g91 Takahashi et al. [TUM Kirk et al. IHMl 
Tavecchio et al. Kino et al. |2"0T)2"|) . However, in 

quasars - but also, in some radio selected BL Lac ob- 
jects, the external radiation fields may be sufficiently 
dense to dominate the Compton cooling of electrons in 
a jet (Madejski et al. I1999|l . Models exploring this hy- 
pothesis are called External-Radiation-Compton (ERC) 
models. Several sources of photons for external radiation 
field have been considered: direct disc radiation (Dermer 
& Schlickeiser ll993[l : diffuse radiation from broad emission 
line (BEL) region (Sikora et n'. 1 99 1 i: infra-red radiation 
from hot dust (Blazejowski et al. l2000fl : or jet synchrotron 
radiation scattered back to the jet by the external gas 
(Ghisellini & Madau lTSSfy . 

In this paper we present a code developed to sim- 
ulate non-thermal flares produced by thin shells propa- 
gating down a conical jet with relativistic speeds. Such 
shells approximate the geometry of relativistic plasma 
enclosed between the forward and reverse shock fronts 
formed by colliding inhomogeneities in a jet (Sikora et 
al. l2Tir)T1 Sikora & Madejski EEEJl- The code is suffi- 



ciently general to treat radiation processes for any ra- 
dial distribution of magnetic and external radiation fields, 
and for any electron/positron injection function. It in- 
cludes both Comptonization of synchrotron radiation and 
Comptonization of external radiation, and takes into ac- 
count adiabatic losses due to 2D conical expansion of the 
shocked plasma sheets. The main limitations of the present 
version of the code are that Comptonization is treated 
self-consistently only within the Thomson limits, and that 
77-pair production is not included. Since there are obser- 
vational indications that in FSRQ the 7-ray spectra have a 
high energy break in the 4- lOGeV band (Pohl et al. 119971 
Sikora et al. 2002), while synchrotron spectra are usually 
steep in the UV band, the above limitations affect only 
marginally our models of quasars. However, they can be 
significant for BL Lac objects, particularly for those bright 
in the TeV band. 

Assumptions and main features of our model, such as 
the equations describing evolution of the energy distri- 
bution of relativistic electrons, radiation processes, rela- 
tivistic aberration and light travel effects are presented 
in Sect. 2. The numerical code is described in Sect. 3. 
Approximate analytical formulas, expressing the model in- 
put parameters as a function of observables, are derived 
in Sect. 4. In Sect. 5 we present results of application of 
our code to model the outburst of 3C279, the best stud- 
ied quasar in the 7-ray band. The work is summarized in 
Sect. 6. 

2. The Model 

In our model a source of the non-thermal radiation is as- 
sumed to be a thin shell, propagating down the conical 
jet with a constant speed, and radiation is produced by 
relativistic electrons/positrons, injected into a shell within 
a given distance range. This picture approximates the in- 
ternal shock scenario, where shocks are formed following 
collisions of inhomogeneities moving with different radial 
velocities (Sikora et al. 119941 Spada et al. l2001[> . Structure 
and dynamics of such internal shocks is in general very 
complex and depends on a number of parameters, such 
as the relative velocity of inhomogeneities, their densi- 
ties, temperatures, geometry, and total masses. The shock 
structure can be double, single (forward or reverse), or can 
initially it can be double, followed by an evolution into 
a single shock (Daigne & Mochkovitch 1199"%! Bicknell & 
Wagner 2002). Regardless, in the case of inhomogeneities 
which are cold prior to their collision and which have com- 
parable masses and comparable rest frame densities, the 
constant speed of the shocked plasma is a good approxi- 
mation. Small thickness of the shocked plasma shell can 
be justified, since internal shocks are at most mildly rela- 
tivistic (Sikora et al. l2UUl[> . 

2.1. Electron Evolution 

Assuming that electron injection function and energy den- 
sities of magnetic field and of external radiation fields are 
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uniform across the shell, one can follow evolution of elec- 
tron energy distribution by solving the kinetic equation 
for the total population of relativistic electrons, despite 
the fact that each element of the conically diverging shell 
has its own rest frame which is different from other frames. 
That equation can be written in the form (Moderski et 
al. HUM: 



dN 1 
dt 1 



d 1 { 7 dt' 



Q, 



(i) 



where Q = (d/V 7 /di')i n j is the electron injection function, 
7 is the random electron Lorentz factor, dj/dt' is the rate 
of the electron energy losses, dt' = dt/T is the proper time, 
and r is the Lorentz factor of the bulk motion of a shell. 
Within the thin shell approximation, the time derivatives 
can be replaced with derivatives over the radial distance 
r. Using a relation dr = fircTdt' , where (3r — vT^-l/f, 
we obtain 



<97V 7 _ 
dr 

where 

d7 
dr 



d_ 



PrcT 



Nr 



d~; 

df 



dr 



rad 



Q 

c/3 r r' 



3 r' 



(2) 



(3) 



The second term on the right-hand side of Eq. © repre- 
sents the adiabatic losses due to two-dimensional (conical) 
expansion of the shell. 

2.2. Radiative Cooling Rates 

There are three processes contributing to the radiative 
energy losses of electrons: 
synchrotron radiation, 



d~; 



<(/37) 2 ; 



s 3m c c 
Comptonization of synchrotron radiation, 



d~; 



ssc 



4(7t 

3to c c 



(5) 



and Comptonization of external radiation, 



d~; 
df 7 



ERG 



3m c c 



where u' B — B' 2 /8ir is the magnetic energy density, u' s is 
the energy density of the synchrotron radiation field, and 
u' ext is the energy density of the external radiation field, 
all as measured in the source comoving frame. B' is the 
strength of the magnetic field, which can be an arbitrary 
function of the radius. 

Energy density of the synchrotron radiation produced 
by a thin shell is 



(f - t' )5L' s ^ (t> - t' )5L' s 



SV 



\>r 2 sn 



2cr 2 Sn- 3 



(7) 



where t' — t' is the time that elapsed from the beginning 
of the shock formation, A' ~ 2c(t' — 1' ) is the width of the 
synchrotron radiation layer at that time, and 



SL' S = 



SL' S ,dv' . 



(8) 



is the synchrotron luminosity produced in a cell of the 
shell, enclosed within a solid angle Sftj <C Oj = Trip 2 
(see Fig. 1). The monochromatic synchrotron luminosity, 
SL' S v , , is computed using Eq. (fTT)f> and integrated over the 
frequency range from u' s min ~ v' ahs , at which the optical 
thickness due to the synchrotron self-absorption is equal 
to unity, to v' s max = (2e/37rm e c)7^ iax B / . 

Energy density of the external radiation field is 



1 



c 

E! 

c 



1 fir 



= - I'dtr = - -^dn 



J 7 cx t(l-/3rcos0 in ) 2 dO 



(9) 



where 6 m is the angle between the trajectory of the in- 
coming photon and the direction of the cell motion. For 
narrow jets this angle can be approximated by the angle 
between the trajectory of the incoming photon and the jet 
axis. 

2.3. Intrinsic Luminosities 
2.3.1. Synchrotron Radiation 

The rate of the synchrotron radiation production per 
cell for a given electron distribution is (Chiaberge & 
Ghisellini lT^jl : 



SL' sy (r) 



(W 7 (r)Fs(i/,7)d7, 



(10) 



(4) where <5iV 7 = iV 7 (£fij/fij), and 



, 3\/3 a T cu' B 2 

^s(^ ,7) = r^x x 



(6) & 



x {K 4/3 (x)K 1/3 ( X ) - -X [K 2 i/3 (x) - Kf /3 ( X )\ j (11) 

is the single electron synchrotron emissivity averaged 
over an isotropic distribution of pitch angles (Crusius 
Schlickeiser 1986), where x = 1/1 / (Sj 2 v B ) and v B = 
eB' j (27rm e c) . The power Fs(p',j) is taken with the low 
energy cutoff at v' ahs . 

Note that because the time scale of the electron/po- 
sitron gyration in the local tangled magnetic fields is 
much shorter than the time scale of radiative energy losses 
(Dermer & Schlickeiser 1993), we can assume that the 
electron momentum distribution is isotropic in the cell co- 
moving frame. It does not imply that the synchrotron radi- 
ation field is isotropic. In general, one must transform the 
synchrotron radiation from neighbor cells into cell frame 
taking into account the change in frequency and direc- 
tion of the incoming photons. However, for a thin shell 



4 



Moderski et al.: Radiation from Blazar Jets 



and small opening angles of the jet the assumption about 
isotropy of the synchrotron radiation field is a reasonable 
approximation. 

2.3.2. Synch rotron-Self-Compton 

Using 5-function approximation, one can find (Chiang & 
Dermer P 



(12) 



XT' ( \ ^ hc(J T /1/2 „ 

SL SSC,u'( r ) - 4 v X 

where 
v[ = Max 



3i/ 



^ In 



and 

vL = Min 



3i/ 3(m c c 2 ) 2 



S,max' 



4 7 2 ' 4/iV 



(13) 



(14) 



where rig ; is the number density of the synchrotron pho- 
tons per frequency. Noting that for isotropic synchrotron 



radiation field 



u' s v , /hu' and 



2r 2 Sn jC 



(15) 



(see Eq. ©), we finally obtain 

sr T l V3(TT ,1/2 

OLaer ,./ = — — — V X 




XT , r-3/2 



dvt. 



(16) 



This is analogous to the Eq. (27) of Chiang & Der- 
mer 



2.3.3. Externa l-Radiation-Compton 

In the co-moving frame, the external diffuse radiation field 
is strongly anisotropic. Due to relativistic aberration, the 
external radiation appears to the jet as mostly incident 
from the forward direction. Within the jet, a photon scat- 
tered by a relativistic electron follows the direction of mo- 
tion of that electron. Because of this, the observer will 
detect preferentially photons that were emitted by the 
electrons which during the scattering process had their 
momentum vector pointing at the observer. In such an 
approximation, the energy of the photon scattered into 
the observer direction is given by (Reynolds [1982 ): 



(17) 



where v' 



rV cxt , hv cxt is the characteristic energy of 



photons in the external radiation field (assumed to have 
a narrow spectrum), and 9' is the angle between the bulk 
motion direction of a given shell segment and direction to 



the observer, as measured in the co-moving frame of the 
bulk motion. For 7 ^> 1 the rate of energy losses of an 
electron with a velocity vector instantaneously oriented 
into direction of the observer is: 



dt> 



ERC 



<7t 
1 

m e c 



4 t (l + cos0') 2 7 2 



(18) 



Note that the angle-averaged value of (d7/d£')ER.c[0'] 
gives Eq. © with /3 ~ 1. Since the "blazar phenomenon" 
implies observing the jet at very small angles, 9 <C tt/2, 
the co-moving observing angle 9' is much smaller than 
7t — 1/r and then (Dermer 1995 ) 



r(l + cos 6>') ~ r(l + /3 r cos6»') = V, 



(19) 



With that approximation the Eqs 1)17(1 and (|18|> can be 

rewritten in the form 



v ~ X>7 v cxt , 
and 



(I7 



dt> 



ERC 



"ext7 

m c c 



(20) 



(21) 



The power per unit solid angle and frequency, emitted by 
a cell via the ERC process into the observer direction is, 
as measured in the cell co-moving frame, 



on' 



-in 



d') 



dt' 



ER.C 



dv' 



(22) 



and after combining with Eqs pUj). J2U, and ©, we ob- 
tain 

^ERcyM cr T (7<5iV 7 )0 



an: 



J (l-/?rcos0 in ) 2 ic Xt dft. 



(23) 



2.4. Observed Spectra and Light Curves 

The monochromatic radiation flux observed at a given in- 
stant is 



TP ft \ - + V" d6L ^ 



where 
d5L v [t,( 



= V 



dl 



,3 wu ^v' V i 



v' = i>/V = v ohs {l + z)/V, 
c(3 r 



cos( 



1 — Pr cos 

cos 9' + p\ 



(t - t ) + r Q . 



1 + [3 r cos9' ' 
t = t ohs /(l + z) , 



(24) 

(25) 
(26) 
(27) 

(28) 
(29) 
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V = T(l + ^ r cos6»') 



1 



T(l -/3 r cos6») 



(30) 



and dh is the luminosity distance. For the Q\ + Q m = 1 
cosmology, 



dx 



o [r> A + a„(l + a;)3]V2 



(31) 



where H is the Hubble constant. 

It should be noted here that because the synchrotron 
and SSC emission is isotropic in the cell co-moving frame, 
the co-moving power per solid angle produced by those 
processes is 



s,ssc 



V ,[Q\ SL S3SCv , 



4tt 



(32) 



where SLsy and SLgscy are given by Eqs l(TU|) and lfT2T) . 

respectively. 

3. Numerical Implementation 

Our numerical method is similar to that used by Chia- 
berge & Ghisellini (1999). Electron evolution equation 
is solved with implicit difference scheme adopted from 
Chang & Cooper H197()[l . First, the uniformly spaced log- 
arithmic energy grid is established: 



Ik = 7o 



In 
7o 



k = 0, . 



(33) 



with energy intervals A^ k = r ) k +i/2 - lk-i/2- Using a 
definition AT? = N 1 ( n f k ,iAr) we can convert Eq. @ into 
a set of linear equations: 



A k N£{ + B k Nl + + C k Nl\\ = Ul 
where the coefficients are 



A k = , 



Cfc 



dr J 



fc-1/2 

Ar 
A 7fc 



Ar 
A7fc 



(34) 



(35) 



fc+l/2 



and the source term is 



ui = ni + 



Ar 



Ql- 



(36) 



The system of equations l|34|) yields a special case of tridi- 
agonal matrix with the upper diagonal equal to zero. It 
is solved numerically by matrix decomposition, forward- 
and back-substitution. In practice we use tridag sub- 
routine from Numerical Recipes (Press et al. I1992|l . This 
method of solving Eq. |(2J) is unstable if the sharp cutoffs 
are present in the electron energy distribution function. In 
our case, because the Eq. © does not contain the acceler- 
ation terms, we avoid this instability by setting the upper 
boundary of the calculation grid to the maximum energy 
of the injected electrons. Another way to avoid sharp cut- 
offs is to introduce very steep power-law tails to mimic 
these cutoffs (Bottcher & Chiang TM?^ . 



There exists an alternative way of solving Eqs (|34|l . If 
one assures that the energy grid is sufficiently wide to 
have N(j n , iAr) = during the entire evolution, then the 
Eq. I|34[) can be solved recursively from the highest to the 
lowest energies using the relation 



C k N, 



k = n — 1, . 



1 



(37) 



For universality, however, we use the full tridag method. 
It allows the consideration of electron distributions with 
high energy tails without performing the calculations for 
unimportant grid points. 

In order to calculate the electron distribution it is nec- 
essary to estimate the electron cooling function ©. This 
requires the knowledge of the energy density of the syn- 
chrotron radiation field which depends on the electron dis- 
tribution itself. In our code this calculation is done itera- 
tively. Initially u' s = is assumed and the electron distri- 
bution is calculated using this value. Then this distribu- 
tion is used to calculate synchrotron radiation density to- 
gether with the frequency of synchrotron self absorption. 
The new value of u' s is used to recalculate the electron 
distribution. This process is repeated until convergence. 
The electron distribution is then used to calculate the syn- 
chrotron luminosity (jTU|) . SSC luminosity ljT2*|) . and ERC 
luminosity (|22|l . For all integrations we use the Romberg's 
method (Press et al. I1992fl . The advantage of this method 
is that it adjusts itself, and thus only a minimal number 
of calculations is performed to achieve desired accuracy. 
Special care must be taken when calculating light curves 
from Eq. 124fl . The radiation reaching the observer at a 
given time is a superposition of radiation emitted at dif- 
ferent radii, and this in turn depends on 9, which is the 
angle between the direction of motion of a given cell and 
the direction to the observer (see Eq. (|27Jl ). 
One must also remember that the choice of the cell size 
puts a constraint on the minimum timescale that can be 
probed with the code. For a given cell size, AO, the min- 
imum timescale for which spatial homogeneity within the 
cell is maintained equals 

Ai~rsin0jA6» (38) 



4. Model Input Parameters 

The input parameters of the model are: 

- the distance of the shock formation from the central 
source, ro; 

- the distance range of shock operation, Ar co n; 

- the bulk Lorentz factor of the shock, T; 

- the energy density of an external diffuse radiation field, 
it cx t and its characteristic frequency ^ ex t; 

- the magnetic field in a shocked plasma, B'; 

- the electron injection function, Q; 

- the angle between the line of sight and the jet axis, V'obs! 

- the half-opening angle of a jet, V'j- 

Initial values of these parameters can be set up by 
using approximate analytical formulas which relate them 
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toward 
observer 



Fig. 1. A schematic picture illustrating the integration 
procedure. A thin shell of electrons enclosed within the 
half-opening angle of the jet, ipj is divided into cells, each 
Srtj thick. Each such cell contains SN^ electrons. The 
observer is located at the angle ipohs from the jet axis. 
Integration is performed by summation over 8 contribu- 
tions from cells located at different radii. 



to the observable quantities. Depending on the available 
observables, a different set of equations must be used. 
Most useful observables used to model the short term 
(1 — 10 days), 7-ray dominated outbursts in FSRQ are: 
the time scale of the outburst/flare, ifl (more specifi- 
cally, its growing part); radiation fluxes at the low-energy 
and high-energy spectral peaks, denoted respectively as 
Fle = v\j&F v ^ and Fhe = ^heFj, he ; location of the spec- 
tral break between the X-ray and 7-ray spectral portions, 
v c ; radiation flux, v+F^, at some v+ > v c ; the spectral 
index of the HE component at v > v c (as measured, e.g., 
by EGRET in its bandpass); and the total radiation flux 
of broad emission lines, Fbel (Celotti et al. H997[l . For the 
sake of simplicity of analytical formulas which will be used 
to derive the initial model parameters as a function of ob- 
servables, both the light travel effects (due to finite size of 
the source) and the Doppler dispersion effects (due to the 
"conical" expansion of a shell) are ignored. In addition, 
we assume that: during the outbursts the production of 
7-rays is dominated by the ERC process; B' oc 1/r; and 
electrons are accelerated in the shock at a constant rate 
and with a single power-law energy distribution, 



Q = K^-v . 

for 7 min = 1 and 



7maX ^ 



<(! + *) 



t^BEL 



(39) 



(40) 



(see Appendix A for justification of such an injection func- 
tion). 

Due to radiative and adiabatic energy losses of elec- 
trons, their energy distribution evolves with time. The 
strong energy dependence of the radiative energy losses 
of electrons causes steepening of the electron energy dis- 
tribution. It takes place at 7 > 7 C , where 7 C is the energy 
at which the time scale of electron energy losses, 



t 



7 



7 



cool 



\dj/dt'\ |<*7M'|rad + \di/dt'U ' 



is equal to the lifetime of a shock 
t&D 



sli 



l + Z 



(41) 



(42) 



During 7-ray dominated outbursts the radiative energy 
losses are dominated by the ERC process, and therefore 



t'cooi = t'sh § ives 

m c c (1 + z)f(k) 
2a T t{iu' t 'D 



(43) 



where f(k) = (1 + 3fc)/(2(l + k)) and k = r /Ar coU . 

At 7 = 7c the power law energy distribution of elec- 
trons, _/V 7 oc 7~ s , changes the slope from s = p at 7 < 7 C 
(slow cooling regime) to s = p + 1 at 7 > 7 C (fast cool- 
ing regime). Such a change of slope results from the fact 
that in the fast cooling regime (i ra d < £sh) the num- 
ber of electrons is saturated by the radiative losses, i.e. 
^7 ~ J 7 Q ^7/M7/^'l ~ QKad- (Note that 7 C is changing 
during the shock progression and that its particular value 
given by Eq. (|43|l is for the instant when the collision is 
complete.) 

The break at 7 C is reflected in the ERC spectral com- 
ponent at 



^c , obs 



7c 2 ^cxt 
l + Z 



(44) 



around which the electromagnetic spectrum changes the 
slope by Aa ~ 0.5. We identify this break as the one 
observed more or less directly in FSRQ between 1 and 
30 MeV. Such identification provides information about 
energy density of external radiation field as measured in 
the shock co- moving frame. Using Eqs l|43|l and l|44|l . we 
obtain 



<xtOf) = 



m e c{l + zfl 2 f{k) 
2cr T tfl(^ c /l'cxt) 1/2 



(45) 



where rf ~ ro + Ar co n ■ Since observations provide directly 
the ratio Fhe/Fle and since 



F HE _ F ERC |d 7 /d<'| ERC [fl'] 
F LE F s ' Id-y/dt'ls 

we can estimate the intensity of the magnetic field, 



u. 
u 



cxt 
/ 

B 



B'{n) 



v 



8ttw' 



Fle 



HE 



(46) 



(47) 



where u' cxt is given by Eq. I|45|l . 
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The short term flares (lasting typically 1-10 days) 
are most likely produced at distances 0.1 — 1.0 pc. At 
such distances, the external photon energy density u' ext 
is dominated by broad emission lines (see, e.g., Sikora et 
al. I2002[l . With this, and assuming spherical geometry of 
the BEL region, we have 



, 4 2 (dL BEL /d\nr)(r) 



(48) 



Provided that (9Lbel/<? In r)(r) is known, this equation, 
combined with Eq. I|45|) and with the formula for the dis- 
tance range of the shock operation, 



Ar co n 



cifl 



i ct fl r 2 v 



1 — Pr cos V'obs 1 + 3 1 + z T 



(49) 



can be used to estimate rf and T. For a broad emis- 
sion line luminosity distribution approximated by (see 
Appendix |B|) 



dL 



BEL 



<91nr 
where 



(r) = 



L 



BEL 



(r/r BEL ) - 5 
(r/r BEL )-°- 5 



for r < r B EL 
for r > tbel 



r BEL - lO^^uv/lO^ergss- 1 ) 7 ^ 



is the distance at which the radial distribution of the BEL 
luminosity has a maximum (Kaspi et al. l2000[) . we obtain 

_ / C?(^c/^BEL)i B EL r BEL if r f < r BEL , . 

r f — 1 2/3/ / U/3 t-2/3 1/3 - t ^ ' y° ' 

Lc r ' (^c/^BEl) 1 ^ BE L r BEL " r f > r BEL 

where i^bel is a characteristic frequency of the broad emis- 
sion lines radiation, and 



(l + z)r t 



V (l + fe)rffl(23/T) 

where c r = cr T (l + z) 1/2 /[67rm e c 3 (l + fc)/(fc)(X>/r)], 
-^bel = 47T(iL-FBEL) and ^uv = 47rd£-Fuv is the luminos- 
ity of the accretion disc. The disc radiation is sometimes 
observed directly (during low states of FSRQ); otherwise 
it can be estimated from Ftjv — -Fbel/£bel, assuming the 
canonical value of the BEL covering factor, £bel = 0.1. 

In order to determine the parameters of the electron 
injection function, we use 7-ray data from the spectral 
band at v 3> v c . There, deeply in the fast cooling regime, 



p = 2a 7 
and 



an 



1 


J 


i 


tot Jl 



Qdj. 



Noting that for 7-ray dominated outbursts, | 7 |tot — 
|7|erCj by using approximate formula for ERC radiation 
production, 



47</Ai^ £Si 0.5[ 7 Ar 7 ]|7|totm e c'r 4 (I?/F) 



one can find that the normalization factor of the electron 
injection function is 

8TTdl[^F Ut } 

where g(a 7 ) — (2a 7 — 1)7^ 2 " T 2 ^ and 




(58) 



Note that for a 7 = 1, g(a 7 ) = 1. 

To complete the set of the model input parameters we 
still need to specify the values of T>/T, k and ipj. Since 
the ERC to synchrotron luminosity ratio is cx (2?/F) 2 
(see Eq. (|4^|) L the view angles V'obs < 1/r are expected 
for strongly 7-ray dominated FSRQ. In such conditions 
T>/Y = 1 ranges from 2 for the observer located on the jet 
axis to 1 for the observer located at -0 o b s = 1/r. 

The value of k is also expected to be of the order 
of unity, as suggested by roughly symmetrical profiles of 

(50) flares (see, e.g., Sikora et al. l2001[) . Regarding the jet open- 
ing angle, we know they are very small (1-3 degrees) on 
kilo-parsec scales, but can be much larger on parsec/sub- 
parsec scales (see, e.g. Lobanov ll998fl . The specific value of 

(51) i/jj can be determined from the value of flux in the soft / mid 
X-ray bands, provided that the latter is dominated by the 
SSC process. Since we do not follow the SSC process in our 
analytical approach, and because other model parameters 
are not very sensitive to i/>j , the numerical simulations can 
be started with any < tpj < 1/r and then corrected in 
subsequent iterations. 

In our method to calculate the model input param- 
eters, we do not take advantage of the ratio of frequen- 
cies, fHE to i/'LE- This ratio often has been used to calcu- 
late the bulk Lorentz factor. If it were to be incorporated 

(53) 

into our method to calculate T, this would give an equa- 
tion for k. However, vF v spectra around their maxima 
are quite flat and therefore, both t^he and i^le are sub- 
ject to very large uncertainties. Furthermore, in FSRQ 
the synchrotron peak is located in the far IR and rarely is 
observed directly. 

An additional difficulty is the fact that the ampli- 
tudes of the 7 -ray flares are usually much larger than the 
amplitudes of the synchrotron flares. This suggests that 
synchrotron component is strongly diluted by the quasi- 
steady radiation produced at larger distances in a jet. We 
can deal with this case in our method by using in the 

(54) Eq. H47JI Ai^ syn instead of F syn , where AF syn is the amount 
by which synchrotron flux increases during the high ampli- 
tude 7 -ray flare. For the instances where the dilution is so 
large that Ai^ syn cannot be determined observationally, as 

(55) an alternative approach, one can use the observed correla- 
tion of the soft /mid X-ray flux with the 7-ray flux. When 
the soft/mid X-ray flux is interpreted as produced by the 
SSC process, one can estimate the magnetic field intensity 
via an iterative process, which relics on the match of the 
model SSC radiation spectrum against that measured in 

(56) the X-ray band. 
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5. Application: blazar 3C 279 

3C 279 is one of the most extensively observed blazars. 
Several times it was a target of multi-wavelength cam- 
paigns. The most fruitful campaign took place in the be- 
ginning of 1996 when the blazar underwent an enormous 
7~ray flare. This event was monitored nearly simultane- 
ously at many frequencies (Wehrle et al. 11998(1 and is an- 
alyzed below using our code. 

5.1. Input parameters 

3C279 is located at a redshift z = 0.538, which for the 
currently favored cosmology (Q,\ ~ 0.7, Q m ~ 0.3 and 
h ~ 0.66) gives the luminosity distance ~ 10 28 cm. 
The data collected during the February 1996 flare (Wehrle 
et al. Hartman et al. l2T)(TT|l give the hv c ~ 20MeV 

and Vi,F Vt ~ 2.75 x 10~ 10 erg s _1 cm -2 , where ~ 
10 23 Hz. The slope of the 7-ray spectrum during the flare 
in the EGRET band was a 7 ~ 0.97 which from Eq. (|5^jl 
gives p ~ 1.9. Although a reliable EGRET spectrum is 
available for this object, the high energy break is not di- 
rectly observed, so we set hu max = 20GeV. Estimate of 
jFbel is very uncertain. Analysis based on relative line in- 
tensities by Celotti et al. |1997§ gives Lbel = 3.4 x 10 44 erg 
s^ 1 . However, during outbursts this luminosity can be 
even larger (Koratkar et al. [1998 ). For our analysis we 
set Lbel = 6.8 x 10 44 crg s _1 . The average photon energy 
of external radiation field is taken to be hv cxt = /i^bel = 
lOeV. Luminosity of the accretion disc, obtained using the 
direct observations of the UV bump during the low state 
of 3C 279 (Pian et al. lT^jl . is L uv ^ 2.7 x 10 45 erg s" 1 . 
Time scale of the flare is t& ~ 1 day (Lawson et al. I1999|l . 

Analysis of the lightcurves of the February 1996 flare 
at different energies (Wehrle et al. I1998[) shows that de- 
spite the very large rise of flux in the 7-ray band, the 
flare was almost undetectable at low energies, especially 
from radio to optical. This suggests that synchrotron com- 
ponent is significantly diluted, and synchrotron radiation 
produced during the flare is hidden by a more steady com- 
ponent. For this reason we initially set F syn / Ferc = 10~ 2 . 
The above set of observables yields (for the assumed 
k = 1 and D/T — 2) the following input parameters: 
< xt (r f ) ~ 0.21erg cm" 3 (<- Eq. g3Jl); B'(r { ) ~ 0.46 
Gauss {*- Eq. (O); r BE L ^ 4.0 x 10 17 cm (<- Eq. JSTJ); 
r t ^ 4.2 x 10 17 cm (<- Eq. fify): T ~ 7.9 (<- Eq. ®); 
K ~ 3.4 x 10 48 s" 1 (<- Eq. (E3J); 7max ~ 3.5 x 10 3 (<- 
Eq.(HOJl). We have also set Vj = 0.12 ~ l/T. 

5.2. Average spectrum 

We first apply our model to the average spectrum of the 
blazar 3C 279 during the February 1996 flare. The result 
is presented in Fig. [3 The fine-tuning of the input param- 
eters requires substantial change only of of the normaliza- 
tion of electron injection function, K. We have achieved 
the best results setting K = 9.5 x 10 49 . This difference can 
easily be understood by examining Eq. 157( 1 . This quan- 
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Fig. 2. The average spectrum of the blazar 3C 279 during 
the February 1996 flare. Data points are from Hartman 
et al. ( 2001). Thick, solid line shows the averaged spec- 
trum of our model (see text for parameters). Thin lines 
represent various components of the spectrum. 

tity is very sensitive to D/T. The calculation of the input 
parameters assumes that the source is point-like and thus 
that the whole radiation is beamed toward the observer. 
Due to the conical structure of the jet the Doppler factor 
decreases to the edge of the jet, and thus more electrons 
are required to produce the observed luminosity. The mi- 
nor correction was applied also for the magnetic field. The 
best result yields B'(rf) = 0.3 Gauss. For the calculation 
we also set ipobs = 0.05, but the averaged spectrum is not 
very sensitive to this value. The averaging was performed 
between r = ro and r = 3ro- 

Large discrepancy between the data and our model in 
the low energy component arises from the fact that the 
synchrotron radiation during the studied flare is substan- 
tially diluted, presumably by radiation produced at larger 
distances in a jet. 

5.3. Light-curves 

Fig. (PJ presents the simulated light-curves together with 
observational data points. Two curves are shown: one is 
the vF v flux at 400 MeV (EGRET range) and the second 
is the integrated 2-20 keV flux (XTE range). For both 
lightcurves the steady components were added. These 
components were estimated from the pre-flare observa- 
tions. In the case of X-ray light curve the quiescent emis- 
sion is F q = 1.23 x 10 _11 erg s _1 cm~ 2 (Lawson et al.Q999) 
and contributes up to 25% of the peak flux. Thus for 
the X-ray band, where the radiation is presumably dom- 
inated by the SSC process, the contribution of the quies- 
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Fig. 3. 7-ray (upper panel - vF v flux at 400 MeV) and Fig. 4. Spectral index a(F v cx i/~ a ) evolution (left pan- 
X-ray (lower panel - integrated 2-10 keV flux) light-curves els) and hardness-intensity diagram (right panels) of the 
of the blazar 3C 279 during the Feb '96 flare. Data points blazar 3C 279 during the February 1996 flare. Upper pan- 
are from Wehrle et al. (1998). Thick, solid line shows the els show the index at 400 MeV while lower panels show 



light-curves of our model. 



the index at 6 keV. 



cent emission is significant. For gamma-rays we estimated 
the steady component from data points before 30 th day of 
1996. Its value is vF UM = 24.7 x 10~ n erg s" 1 cm" 2 and 
its contribution to the peak flux is only 12%. 

The agreement between the data and our model is 
remarkable taking into account the relative simplicity of 
the model and complexity of the real phenomena as sug- 
gested by, e.g., the gamma-ray light-curve. The discrep- 
ancy at the rising part of the X-ray light curve sug- 
gests an existence of the soft X-ray precursors caused by 
Comptonization of the external radiation by cold electrons 
(Sikora & Madeiski lMi^l . 

5.4. Spectral evolution 

In order to demonstrate spectral changes during the flare, 
in Fig. 2] we present the evolution of the spectral index 
a (defined by F v cx v~ a ) at two energies: 400 MeV and 
6 keV. Two plots are presented for each energy: one is 
the power law spectral index a as a function of time, the 
second is the so called hardness-intensity diagram showing 
the spectral index as a function of the flux. 

Unfortunately, the comparison with observations in 
this case is difficult. There are no observations of the 
spectral index variation during the flare by EGRET. The 
reported slope of the 7-ray spectrum for the event was 
a-y — 0.97 which is in excellent agreement with the mid- 
dle part of the curve in upper left corner of the Fig. 0] 
For the XTE observations, a x ranges from 0.68 to 1.26, 
but the errors are large, and there is a possible evidence 



of systematic error of 0.1 (Lawson et al. I1999[) . Even with 
the errors taken into account, the observations seem to 
disagree with results presented in lower left panel of the 
Fig. However, we note that in this range, the influence 
of the quiescent component is significant, and the proper 
comparison would require the analysis which is beyond 
the scope of this paper. The steady component may also 
influence the spectrum at 7-ray energies and cause the dis- 
crepancy for observations outside of the maximum of the 
flare. 

5.5. Evolution of Electron Distribution 

In Fig. we present the time evolution of energy distribu- 
tion of electrons during the flare. The evolution is followed 
from r = ro till the electrons reach the distance 3 ro. At 
the very beginning of the flare, the number of all elec- 
trons increases, then the number of the most relativistic 
electrons starts to saturate being balanced by the radia- 
tive energy losses. After the injection of electrons stops, 
the high energy tail of the electron distribution decays 
very rapidly, while only small changes can be noticed at 
lowest energies. 

The characteristic feature of this evolving electron 
spectrum is the spectral break, which moves from higher 
to lower values as evolution proceeds. At a distance r = 2 
ro, where injection stops and the 7 2 -/V 7 reaches maximum, 
the approximate location of the break is 7 ~ 65. The value 
calculated from Eq. Ij43(l is somewhat larger {p[ c ~ 105) 
and this is because the approximate analytical formula 
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Fig. 5. Evolution of the electron energy distribution dur- 
ing the February 1996 flare of the blazar 3C 279. Step 
between the curves is Ar = 0.2ro. 

for 7 C does not include synchrotron and SSC energy losses. 
At energies lower than 7 C (the slow cooling regime), the 
power-law energy distribution of electrons N 7 ~ 7~ s has 
an index s — p, while for 7 3> 7c (the fast cooling regime) 
the slope of the electron distribution (due to efficient ra- 
diative cooling) is steeper by unity i.e., s = p + 1. The 
latter effect is mostly obscured by the effect of the break 

^t 7max- 

6. Summary 

In this article, we have presented our numerical code 
BLAZAR which simulates light-curves and spectra of 
blazars during flares. In the code, the structure of the 
source responsible for the production of flares is approx- 
imated by thin uniform shells, propagating with con- 
stant speed down the conical jet. The shells are filled 
with relativistic electrons/positrons, which are injected 
within a given distance range at a constant rate. Evolution 
of the electron energy distribution is treated by the ki- 
netic/continuity equation, given by Eq. Q. In that equa- 
tion, the electron injection function Q is separated from 
the cooling term. This simplification can be justified pro- 
vided time scale of acceleration of electrons is much 
shorter than time scale of their cooling and time scale of 
the shock operation. Since in FSRQ, which are main tar- 
gets of the present version of the code, the cooling break 
is located at v c -C i/ ma x, the above assumption is satisfied 
for all electrons except for those with 7 ~ 7max- However, 
this condition may not be satisfied for the X-ray selected 
BL Lac objects, including TeV blazars. In these objects 
time scale of acceleration can well be comparable to the 



lifetime of the source, and then Eq. must be modified 
in the way presented by Kirk et al. I|1994f) (see also Kirk 
et al. rHT9 8 and Bottcher & Chiang 2002 for application to 
blazars). 

Very short gyration time scale (as compared with 
the radiative energy losses) plus very likely presence of 
chaotic magnetic fields (because of MHD turbulence in- 
duced around the shock fronts) justify the assumption 
regarding the isotropic distribution of the electron mo- 
menta, and therefore regarding the isotropy of the syn- 
chrotron and SSC emission in the source co-moving frame. 
However, the isotropy approximation doesn't apply to 
Comptonization of external radiation. This is because due 
to relativistic aberration, the external diffuse radiation is 
seen by the source as strongly beamed from the front. 
Because of this, the scattered radiation is anisotropic - it is 
beamed into the source propagation direction already even 
in the source co-moving frame. This anisotropy, originally 
pointed out by Dermer fTjJ55$, is self-consistently treated 
in our code. Radiative effects of the ERC anisotropy, of- 
ten ignored in other models, are particularly strong for 

V'obs > 1/r. 

Another advantage of our code is that it includes the 
adiabatic energy losses in the kinetic equation for elec- 
trons. In our model they are related to the "conical" ex- 
pansion of a shell. Inclusion of adiabatic losses is very im- 
portant in calculating properly the location of the cooling 
break v c . Those losses also strongly affect the light curves 
of flares observed at v < v c (Sikora ct al. 2001). 

As an example of application of our code, we presented 
in §5 the results of our modeling of the short term out- 
burst observed in 3C 279 in February 1996 (Wehrle et 
al. 1998). Our results demonstrate very strong dilution of 
the synchrotron radiation component in this object. This 
effect, often ignored completely by other models, can be 
quite common, as suggested by much smaller amplitudes 
of optical flares than of 7-ray flares (Ulrich et al. ^997). 
For such objects, detailed fitting of the synchrotron spec- 
trum in terms of the homogeneous/one-component model 
is inappropriate and can lead to very large errors in the 
model parameters. Our method, albeit approximate, takes 
into account the dilution effect. 

It is worth noting that the code "BLAZAR" was al- 
ready used in several previous papers. It was applied 
to the study the dependence of amplitude of flares in 
blazars on the frequency, on the ratio r$ / Ar co u, and on the 
time profile of an electron injection function Q(t') (Sikora 
et al. I2001fl . The code was also used to study the relative 
role of hot dust and broad emission line region as a source 
of the external seed photons for the inverse Compton pro- 
cess in a jet (Blazejowski et al.) and was applied to demon- 
strate the possible unification between the MeV-blazars 
and GeV-blazars (Sikora et al. 2002|) . Finally, after some 
small modifications allowing to model the radiation from 
the shells propagating with a variable bulk Lorentz factor 
and supplementing the code by relevant dynamical equa- 
tions, it was used to study the effect of lateral expansion 
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of the collimated GRBs, possibly resulting in the observed 
breaks in their light-curves (Moderski et al.|5000). 

In the present form, the code doesn't include the Klein- 
Nishina effects nor the 7-ray absorption and e + e~ pair 
creation. In FSRQ the Klein-Nishina effects become im- 
portant at 7 > a few x 10 4 /T, while the absorption - at 
hv > 20GeV, both due to interactions with external BEL 
photons. The high energy break at hv ~ 4 — lOGeV, sug- 
gested by EGRET observations of FSRQ during the 7-ray 
flares(Pohl et al. 1997), and the steep UV spectra seem 
to indicate that energy distribution of electrons/positrons 
has a break or cutoff around 7 ~ a few x 10 4 /r (see 
Appendix A). Hence, the KN and 77 pair production 
effects are expected to be at most marginal. The above 
limitations, particularly the first one, are expected to be 
relevant only for BL Lac objects, in particular for those 
with strong TeV emission. In the case of our model of 
the 3C 279 flare T^ max hv^i J /m e c 2 ~ 0.5 and hv maK = 
20GeV, so the effects are expected to be negligible. 
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Appendix A: Electron injection function 

In most FSRQ, the high energy spectra of flares can be 
recovered assuming a single power-law injection function. 
Due to cooling effect, the electron energy distribution 
evolves into the broken power-law, which is reflected in 
the spectrum of the ERC component, with the break lo- 
cated in the 1-30 MeV range (see §5). The low energy tails 
observed in several FSRQ do not show any change of the 
slope down to the lowest X-ray energies, and this seems 
to indicate that the power-law electron injection func- 
tion extends down to 7 min <few (Tavecchio et al. 119981 
Sikora et al. 2002). This is independently supported by 
the detection of circular polarization in radio cores of some 
quasars (Homan et al. 12001(1 . provided such polarization 
is generated by the Faraday conversion process (Wardle 
et al.GIIIIHI Ruszkowski & Begelman lM)^ . 

It is less clear how electron injection function behaves 
at the highest energies and what is the maximum en- 
ergy of injected electrons. EGRET data of FSRQ during 
flares suggest existence of a high-energy break at 4 — 10 
GeV (Pohl et al. I1997[) . This high-energy break cannot 
be caused by the Klein- Nishina effect, because in the fast 
cooling regime, the luminosity of the Compton component 
is determined by the electron injection function and reduc- 
tion of the Compton scattering cross-section in the Klein- 
Nishina regime can be substantially compensated by the 
respective increase of the electron densities. In particular, 
for electron injection function oc 7~ 2 the Klein-Nishina 
portion of the Compton spectrum has exactly the same 
spectral index, a = 1, as in the Thomson regime (Zdziarski 
& Krolik 1993). The high-energy break also cannot be pro- 
duced by absorption of 7-rays by external photons. This is 
because the radiative environment is transparent for pho- 
tons detected within the EGRET band (it becomes opaque 
only for hv > 20 GeV, due to absorption of such photons 
by optical/UV broad emission lines). Hence, the break in 
4—10 GeV band, if real, must be related to the break in 
the electron injection function, at 7 ~ a few x 10 4 /r. 

Is the injection function extending much above that 
break, but with a steeper slope as suggested by observa- 
tions of the so called MeV-blazars (Sikora et al. PJ002), 
or, there is a cutoff? One can try to answer that ques- 
tion by studying the spectra in the UV band, where the 
high energy tails of the synchrotron component dominate. 
However, the task is not easy, because synchrotron radi- 
ation in this band can be affected by several factors: by 
extinction in the host galaxy; by flattening of energy dis- 
tribution of highest energy electrons due to Compton scat- 
terings in the KN regime (Dermer & Atovan [2002l) and by 



77 pair production; and, as in the case of 3C 279, by dilu- 
tion of the synchrotron component by radiation produced 
at other locations in a jet. Hence, the possible extension 
of the injection function up to energies 7 3> a few x 10 4 /r 
can be verified observationally in FSRQ only after the 
launch of the GLAST satellite in 2006. 

Fortunately, the basic features of the current model 
of FSRQ as adopted here are not very sensitive to the 
presence and details of steep high energy tails of the in- 
jection function, and, in order to reproduce the spectra 
at the level of detail required by the current data, it is 
sufficient to use a single power-law injection function as a 
first approximation. 

Appendix B: Broad emission line region 

Reverberation mappings indicate that the BEL region 
in AGNs is geometrically thick, with r max /r m i n > 30 
( Peterson IM)3) . The 3D structure of the BEL region is 
still not fully understood, and in our code is assumed to be 
spherical. We approximate the radial distribution of BEL 
luminosity by two power-law functions which join at tbel 
where the luminosity has a peak, i.e. 
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and qi > and q2 > 0. In the case of a smooth transition 
from the BEL region to the narrow emission lines region, 
one can expect qi ~ 0.5. The value of qi is very uncertain 
and, in order to minimize number of parameters, we adopt 
qi = q 2 = 0.5. 



